class: center, middle, inverse, title-slide # Espace et Temps avec R ## 1 - Manipulation de données spatiales ### Robin Cura & Lise Vaudor ### 16/10/2018
École Thématique GeoViz 2018 --- ## Sommaire .pull-left[ Charger des données spatiales - [Les Simples Features et R](#objets-sf) - [Lecture de fichiers spatiaux](#lecture-sf) - [Conversion de tableaux en sf](#conversion-spatiale) Manipuler des données spatiales - [Éditer des données](#modifier-sf) - [Systèmes de projections](#projections-sf) - [Agrégations spatiales](#agregations-sf) - [Jointures spatiales](#jointures-spatiales) - [Opérations géométriques](#operations-spatiales) - [Quelques exemples d'analyse spatiale avec R](#analyse-spatiale) ] .pull-right[ Exploration visuelle de données spatiales - [Afficher un objet avec plot](#plot-sf) - [Exploration rapide avec mapview](#plot-mapview) - [Exploration thématique avec mapview](#carto-mapview) Cartographier avec R - [Cartographier avec ggplot2](#carto-ggplot2) - [Cartographie statique avec cartography](#carto-cartography) - [Cartographie dynamique avec leaflet](#carto-leaflet) ] --- ## Les données spatiales dans R Deux formats de données spatiales coexistent : - `sp`, le format le plus ancien et répendu - `sf`, un format plus récent, puissant, mais pas encore universel .pull-left[ .smaller[ ```r library(rgdal) library(sp) iris_Paris_sp <- readOGR("data/IRIS_Paris_L93.shp", stringsAsFactors = FALSE) ``` ``` ## OGR data source with driver: ESRI Shapefile ## Source: "/data/user/c/rcura/geoviz2018/data/IRIS_Paris_L93.shp", layer: "IRIS_Paris_L93" ## with 987 features ## It has 5 fields ``` ]  ] .pull-right[ .smaller[ ```r library(sf) iris_Paris_sf <- st_read("data/IRIS_Paris_L93.shp", stringsAsFactors = FALSE) ``` ``` ## Reading layer `IRIS_Paris_L93' from data source `/data/user/c/rcura/geoviz2018/data/IRIS_Paris_L93.shp' using driver `ESRI Shapefile' ## Simple feature collection with 987 features and 5 fields ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: 645053.4 ymin: 6857483 xmax: 657162.5 ymax: 6867081 ## epsg (SRID): NA ## proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs ``` ```r str(iris_Paris_sf) ``` ``` ## Classes 'sf' and 'data.frame': 987 obs. of 6 variables: ## $ Dept : chr "75" "75" "75" "75" ... ## $ DepCom : chr "75101" "75101" "75101" "75101" ... ## $ Nom_Com : chr "PARIS 1ER" "PARIS 1ER" "PARIS 1ER" "PARIS 1ER" ... ## $ Iris : chr "0101" "0102" "0103" "0104" ... ## $ DComIris: chr "751010101" "751010102" "751010103" "751010104" ... ## $ geometry:sfc_MULTIPOLYGON of length 987; first list element: List of 1 ## ..$ :List of 1 ## .. ..$ : num [1:21, 1:2] 652145 652121 652096 652078 651998 ... ## ..- attr(*, "class")= chr "XY" "MULTIPOLYGON" "sfg" ## - attr(*, "sf_column")= chr "geometry" ## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA ## ..- attr(*, "names")= chr "Dept" "DepCom" "Nom_Com" "Iris" ... ``` ] ] --- ## Les données spatiales dans R .pull-left[ .smaller[ ```r plot(iris_Paris_sp) ``` <!-- --> ] ] .pull-right[ .smaller[ ```r plot(iris_Paris_sf) ``` <!-- --> ] ] --- ## Les données spatiales dans R On va se contenter de mobiliser le *package* `sf`, qui permet une manipulation d'objets simples, semblables à des `data.frame`/`tibbles`. **N.B.** : Pour une très grande majorité d'opérations, `sf` suffit. Il peut toutefois être nécessaire de passer de `sf` à `sp` dès lors qu'on souhaite utiliser des fonctions spatiales avancées, par exemple pour mobiliser des méthodes d'analyse spatiale. --- name: objets-sf ## `sf` : Les *Simples Features* .small[ .center[  R. Lovelace, J. Nowosad & J. Muenchow (2018), Geocomputation with R - https://geocompr.robinlovelace.net/ ] ] --- ## `sf` : Structure d'un objet `sf` .small[ ```r iris_Paris_sf ``` ``` ## Simple feature collection with 987 features and 5 fields ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: 645053.4 ymin: 6857483 xmax: 657162.5 ymax: 6867081 ## epsg (SRID): NA ## proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs ## First 10 features: ## Dept DepCom Nom_Com Iris DComIris geometry ## 1 75 75101 PARIS 1ER 0101 751010101 MULTIPOLYGON (((652144.5 68... ## 2 75 75101 PARIS 1ER 0102 751010102 MULTIPOLYGON (((651758.7 68... ## 3 75 75101 PARIS 1ER 0103 751010103 MULTIPOLYGON (((651760.1 68... ## 4 75 75101 PARIS 1ER 0104 751010104 MULTIPOLYGON (((651564.6 68... ## 5 75 75101 PARIS 1ER 0105 751010105 MULTIPOLYGON (((650420.6 68... ## 6 75 75101 PARIS 1ER 0199 751010199 MULTIPOLYGON (((652068.4 68... ## 7 75 75101 PARIS 1ER 0201 751010201 MULTIPOLYGON (((652342.9 68... ## 8 75 75101 PARIS 1ER 0202 751010202 MULTIPOLYGON (((652091 6862... ## 9 75 75101 PARIS 1ER 0203 751010203 MULTIPOLYGON (((652041.8 68... ## 10 75 75101 PARIS 1ER 0204 751010204 MULTIPOLYGON (((652228.7 68... ``` ] - `\(\rightarrow\)` Un `data.frame`, doté d'une colonne de `geometry` et d'un **système de coordonnées/projection**. --- name: lecture-sf ## Lecture de fichiers géographiques #### Toutes les fonctions de `sf` s'inspirent de la syntaxe des fonctions `PostGIS` : `st_OPERATION` : - Lecture d'un shapefile .small[ ```r irisIDF <- st_read(dsn = "data/IRIS_Paris_L93.shp", stringsAsFactors = FALSE) # Comme pour la lecture des data.frame ``` ``` ## Reading layer `IRIS_Paris_L93' from data source `/data/user/c/rcura/geoviz2018/data/IRIS_Paris_L93.shp' using driver `ESRI Shapefile' ## Simple feature collection with 987 features and 5 fields ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: 645053.4 ymin: 6857483 xmax: 657162.5 ymax: 6867081 ## epsg (SRID): NA ## proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs ``` ] --- name: projections-sf ## Projections / transformations #### Comme dans tout SIG, le système de projection (`st_crs`) est récuperé à la lecture d'un objet : .small[ ```r st_crs(irisIDF) ``` ``` ## Coordinate Reference System: ## No EPSG code ## proj4string: "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs" ``` ] Ici, le "code de projection" est bien reconnu, mais le code EPSG n'est pas renseigné. On va le spécifier pour clarifier (même si ce n'est pas utile) : .small[ ```r irisIDF <- irisIDF %>% st_set_crs(2154) # SRID/EPSG du Lambert 93 head(irisIDF) ``` ``` ## Simple feature collection with 6 features and 5 fields ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: 650184.6 ymin: 6861751 xmax: 652176.1 ymax: 6863139 ## epsg (SRID): 2154 ## proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ## Dept DepCom Nom_Com Iris DComIris geometry ## 1 75 75101 PARIS 1ER 0101 751010101 MULTIPOLYGON (((652144.5 68... ## 2 75 75101 PARIS 1ER 0102 751010102 MULTIPOLYGON (((651758.7 68... ## 3 75 75101 PARIS 1ER 0103 751010103 MULTIPOLYGON (((651760.1 68... ## 4 75 75101 PARIS 1ER 0104 751010104 MULTIPOLYGON (((651564.6 68... ## 5 75 75101 PARIS 1ER 0105 751010105 MULTIPOLYGON (((650420.6 68... ## 6 75 75101 PARIS 1ER 0199 751010199 MULTIPOLYGON (((652068.4 68... ``` ] --- ## Projections / transformations #### On peut re-projeter un objet `sf` avec la fonction `st_transform(SRID/CRS)` : .pull-left[ .small[ ```r irisIDF %>% select(geometry) %>% # Lambert 93 plot() ``` <!-- --> ] ] .pull-right[ .small[ ```r GallPeters <- "+proj=cea +lon_0=0 +lat_ts=45 \ +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs" irisIDF %>% select(geometry) %>% * st_transform(crs = GallPeters) %>% plot() ``` <!-- --> ] ] --- name: plot-sf ## Visualisation exploratoire : `plot()` On peut utiliser la fonction de base `plot()`, qui affiche alors l'ensemble des attributs du jeu de données : .small[ ```r plot(irisIDF) ``` <!-- --> ] --- ## Visualisation exploratoire : `plot()` Pour "cartographier" le contenu d'une seule variable : .small[ ```r plot(irisIDF["Nom_Com"]) ``` <!-- --> ] --- name: plot-mapview ## Visualisation exploratoire : `mapview` Le *package* `mapview` permet de mener rapidement une exploration des données spatiales et attributaires : .small[ ```r library(mapview) # install.packages("mapview") mapview(irisIDF) ```
] --- ## Visualisation exploratoire : `mapview` Le *package* `mapview` permet de mener rapidement une exploration des données spatiales et attributaires. On peut aussi spécifier un attribut à observer avec l'instruction `zcol` : .small[ ```r mapview(irisIDF, zcol = "Nom_Com") ```
] --- name: modifier-sf ## Opérations attributaires Comme l'objet `sf` est fondamentalement un `data.frame`, on peut lui appliquer les manipulations vues hier : filtrage (`filter`), réductions de variables (`select`), renommage de variables (`rename`) et création/modification de variables (`mutate`) : .small[ ```r irisIDF %>% * select(DepCom, DComIris) %>% # La sélection de la colonne `geometry` est implicite * mutate(ARRONDISSEMENT = as.numeric(DepCom) - 100) %>% * filter(ARRONDISSEMENT >= 75001, ARRONDISSEMENT <= 75007) %>% mapview(zcol = "ARRONDISSEMENT") ```
] --- name: agregations-sf ## Agrégations Comme avec un `tibble`, on peut réaliser des opérations d'agrégation. La géométrie est alors, elle-aussi, agrégée en conséquence : .small[ ```r arrondissements_Paris <- irisIDF %>% group_by(DepCom) %>% summarise() plot(arrondissements_Paris) ``` <!-- --> ] --- ## Jointures attributaires Comme pour tous les `data.frame`, on peut réaliser des jointures attributaires : .smaller[ ```r df_dmr <- readRDS("dans_ma_rue_clean.RDS") joinData <- df_dmr %>% group_by(CODE_POSTAL, ANNEE_DECLARATION) %>% summarise(NbIncidents = n()) %>% ungroup() donnees_incidents <- arrondissements_Paris %>% mutate(CODGEO = as.character(as.numeric(DepCom) - 100)) %>% * left_join(joinData, by = c("CODGEO" = "CODE_POSTAL")) incidents2015 <- filter(donnees_incidents, ANNEE_DECLARATION == 2015) incidents2017 <- filter(donnees_incidents, ANNEE_DECLARATION == 2017) commonBreaks <- joinData %>% filter(ANNEE_DECLARATION %in% c(2015, 2017)) %>% pull(NbIncidents) %>% quantile(probs = seq(from = 0, to = 1, by = 0.1)) ``` .pull-left[ ```r plot(incidents2015[,"NbIncidents"], main = "Incidents 2015", breaks = commonBreaks) ``` <!-- --> ] .pull-right[ ```r plot(incidents2017[,"NbIncidents"], main = "Incidents 2017", breaks = commonBreaks) ``` <!-- --> ] ] --- name: conversion-spatiale ## Conversion en `sf` et jointures spatiales Les objets `sf` sont des objets spatiaux, on peut donc aussi effectuer des jointures spatiales (entre deux objets `sf`) : - Conversion du jeu de données "Dans ma Rue" en objet spatial avec la fonction `st_as_sf` : .small[ ```r dmr_spatial <- df_dmr %>% filter(ANNEE_DECLARATION %in% c(2015, 2017)) %>% * st_as_sf(coords = c("Long", "Lat"), crs = 4326) %>% st_transform(2154) # Pour une opération géométrique, les objets doivent avoir le même SRID ``` ] --- name: jointures-spatiales ## Jointures spatiales - Jointure spatiale pour récuperer les IRIS contenant les points : .small[ ```r dmr_augmente <- dmr_spatial %>% * st_join(irisIDF %>% select(DComIris), join = st_within) ``` .pull-left[ ```r dmr_augmente %>% group_by(DComIris, ANNEE_DECLARATION) %>% summarise(NbIncidents = n()) %>% ungroup() %>% head() ``` ] ] .smaller[ .pull-right[ ``` ## Simple feature collection with 6 features and 3 fields ## geometry type: MULTIPOINT ## dimension: XY ## bbox: xmin: 651589.7 ymin: 6861822 xmax: 652155.5 ymax: 6862494 ## epsg (SRID): 2154 ## proj4string: +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ## # A tibble: 6 x 4 ## DComIris ANNEE_DECLARATI… NbIncidents geometry ## <chr> <int> <int> <MULTIPOINT [m]> ## 1 751010101 2015 52 (651770.5 6862257, 651770.8 6862… ## 2 751010101 2017 61 (651802.6 6862245, 651834.4 6862… ## 3 751010102 2015 6 (651690 6862090, 651709.4 686207… ## 4 751010102 2017 10 (651674.3 6862086, 651679.7 6862… ## 5 751010103 2015 14 (651589.7 6862336, 651599.1 6862… ## 6 751010103 2017 15 (651620.5 6862408, 651637.9 6862… ``` ] ] --- ## Opérations spatiales name: operations-spatiales `sf` permet de réaliser la quasi-totalité des opérations attendues d'un SIG, avec des fonctions dédiées reprenant le vocabulaire classique (inspiré par `PostGIS`) : `st_agr`, `st_area`, `st_bbox`, `st_bind_cols`, `st_boundary`, `st_buffer`, `st_cast`, `st_centroid`, `st_combine`, `st_contains`, `st_contains_properly`, `st_convex_hull`, `st_covered_by`, `st_covers`, `st_crop`, `st_crosses`, `st_difference`, `st_disjoint`, `st_distance`, `st_equals`, `st_equals_exact`, `st_graticule`, `st_interpolate_aw`, `st_intersection`, `st_intersects`, `st_is_simple`, `st_is_valid`, `st_is_within_distance`, `st_jitter`, `st_join`, `st_layers`, `st_length`, `st_line_merge`, `st_line_sample`, `st_make_grid`, `st_overlaps`, `st_point_on_surface`, `st_polygonize`, `st_relate`, `st_sample`, `st_segmentize`, `st_simplify`, `st_snap`, `st_sym_difference`, `st_touches`, `st_transform`, `st_triangulate`, `st_union`, `st_viewport`, `st_voronoi`, `st_within`, `st_wrap_dateline`, `st_write`, `st_write_db`, `st_zm`, `write_sf` Voir les exemples illustrées dans les *vignettes* du *package* `sf` : [Manipulating Simple Feature Geometries](https://cran.r-project.org/web/packages/sf/vignettes/sf3.html) et [Manipulating Simple Features](https://cran.r-project.org/web/packages/sf/vignettes/sf4.html) --- ## Opérations spatiales : un exemple de traitement -> On peut chercher à créer une carte des densités d'incidents déclarés par IRIS, potentiellement selon les années : .small[ ```r nbIncidentsIris <- dmr_spatial %>% st_join(irisIDF %>% select(DComIris), join = st_within) %>% st_set_geometry(NULL) %>% # On n'a plus besoin que ce soit un objet spatial group_by(DComIris, ANNEE_DECLARATION) %>% * summarise(NbIncidents = n()) %>% ungroup() irisIDF_Incidents <- irisIDF %>% * mutate(surface = st_area(.)) %>% left_join(nbIncidentsIris, by = "DComIris") densite_incidents <- irisIDF_Incidents %>% * mutate(densiteIncidents = NbIncidents / (surface / 1E6)) %>% # Pour avoir une densité / km² select(DepCom, DComIris, ANNEE_DECLARATION, NbIncidents, densiteIncidents) # Les densités sont exprimées avec des "unités", qui complexifient le traitement. # On les enlève donc : densite_incidents <- densite_incidents %>% mutate(densiteIncidents = units::drop_units(densiteIncidents)) ``` ] --- name: carto-mapview ## Cartographie exploratoire : `mapview` .small[ ```r densite2015 <- filter(densite_incidents, ANNEE_DECLARATION == 2015) densite2017 <- filter(densite_incidents, ANNEE_DECLARATION == 2017) map2015 <- mapview(densite2015, zcol = "densiteIncidents", legend = TRUE) map2017 <- mapview(densite2017, zcol = "densiteIncidents", legend = TRUE) latticeView(map2013, map2017, ncol = 2) ``` ] --- ## Cartographie exploratoire : `mapview` .small[
] --- name: carto-ggplot2 ## Cartographie statique : `ggplot2` On peut bien sûr afficher ces éléments avec le *package* `ggplot2`, en faisant appel à la géométrie dédiée aux objets de type `sf` : `geom_sf` .small[ ```r ggplot(densite_incidents) + geom_sf(data = densite_incidents, aes(fill = densiteIncidents), lwd = 0) + * # On est obligé de remettre la couche sf en data facet_wrap(~ANNEE_DECLARATION) + scale_fill_viridis_c(name = "Incidents / Km²") + guides(fill = guide_colourbar(title.position = "top")) + * coord_sf() + # On s'assure que les coordonnées du graphique respectent le CRS theme_minimal() + theme(legend.position = "bottom", legend.key.width = unit(2,"cm")) ``` <!-- --> ] --- ## Cartographie statique : `ggplot2` - Certains *packages* (`ggsn`, `ggmap`) proposent des foncitonnalités permettant d'améliorer le rendu cartographique (échelles graphiques, flèches d'orientations, fonds de cartes...) ; - mais `ggplot` est avant tout un outil de visualisation **générique** : #### - pour créer des cartographies de qualité, mieux vaut se tourner vers des *packages* dédiés : `tmap` (qu'on ne verra pas ici) et `cartography` --- name: carto-cartography ## Cartographie statique : `cartography` #### Exemple d'une carte en symboles proportionnels .small[ ```r library(cartography) # On commence par afficher un fond de carte par(mar=c(0,0,0,0)) plot(st_geometry(arrondissements_Paris), col = "#1F4257", border = "#B6C1C8", bg = "#708694", lwd = 2) # On ajoute des symboles proportionnels pour les incidents propSymbolsLayer(x = densite_incidents %>% filter(ANNEE_DECLARATION == 2017), var = "NbIncidents", col = "#F97B64", inches = 0.1, legend.title.txt = "Nombre de déclarations par IRIS en 2017") # On ajoute les labels des arrondissements labelLayer(x = arrondissements_Paris %>% mutate(DepCom = as.numeric(DepCom) - 100), txt = "DepCom", col = "#B6C1C8", bg = "#1F4257", halo = TRUE, overlap = FALSE, show.lines = FALSE) ``` ] --- ## Cartographie statique : `cartography` #### Exemple d'une carte en symboles proportionnels .center[ .small[ <!-- --> ] ] --- name: carto-leaflet ## Cartographie dynamique : `leaflet` #### - Pour créer des cartes "dynamiques", en plus de `mapview` (qui sert surtout à l'exploration de données), on peut utiliser le *package* `leaflet`, qui génère une carte en HTML+JavaScript que l'on pourra placer sur une page web quelconque : .small[ ```r library(leaflet) densite_incidents2017 <- densite_incidents %>% filter(ANNEE_DECLARATION == 2017) %>% st_transform(4326) # Leaflet requiert une couche en WGS84 seuils <- quantile(densite_incidents2017$densiteIncidents, probs = seq(from = 0, to = 1, by = 0.2)) colorPalette <- colorBin("YlOrRd", domain = densite_incidents2017$densiteIncidents, bins = seuils) infosPopup <- sprintf("<strong>%s</strong><br/>%s incidents (%.1f / km²)", densite_incidents2017$DComIris, densite_incidents2017$NbIncidents, densite_incidents2017$densiteIncidents) %>% lapply(htmltools::HTML) leaflet(data = densite_incidents2017) %>% addProviderTiles(providers$CartoDB.DarkMatterNoLabels) %>% addPolygons( fillColor = ~colorPalette(densiteIncidents), fillOpacity = 0.7, color = "white", weight = .5, opacity = 1, dashArray = "3", label = infosPopup) %>% addLegend(pal = colorPalette, values = ~densiteIncidents, opacity = 0.7, title = "Densités d'incidents<br/>[incident/km²]", position = "topright") ``` ] --- ## Cartographie dynamique : `leaflet` #### - Pour créer des cartes "dynamiques", en plus de `mapview` (qui sert surtout à l'exploration de données), on peut utiliser le *package* `leaflet`, qui génère une carte en HTML+JavaScript que l'on pourra placer sur une page web quelconque : .small[
] --- name: analyse-spatiale ## Un peu d'analyse spatiale ?